Solving real time evolution problems by constructing excitation operators 
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In this paper we study the time evolution of an observable in the interacting fermion systems 
driven out of equilibrium. We present a method for solving the Heisenberg equations of motion by 
constructing excitation operators which are defined as the operators A satisfying [H, A] — XA. It 
is demonstrated how an excitation operator and its excitation energy A can be calculated. By an 
appropriate supposition of the form of A we turn the problem into the one of diagonalizing a series 
of matrices whose dimension depends linearly on the size of the system. We perform this method 
to calculate the evolution of the creation operator in a toy model Hamiltonian which is inspired by 
the Hubbard model and the nonequilibrium current through the single impurity Anderson model. 
This method is beyond the traditional perturbation theory in Keldysh-Green's function formalism, 
because the excitation energy A is modified by the interaction and it will appear in the exponent in 
the function of time. 

PACS numbers: 02.70.-c, 72.15.Qm 
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I. INTRODUCTION 

The interacting fermion systems have a central position 
in modern condensed matter physics [l[ . Their properties 
have been intensively studied both in theory and experi- 
ment. Good agreements between the theoretical predic- 
tions and the experimental results have been obtained 
when these systems are in thermal equilibrium. Exam- 
ples include the Coulomb blockade effect and the Kondo 
effect in quantum dots 0-Q > and the charge density wave 
in the low-dimensional materials, e.g. the carbon nan- 
otubes and the edge states in fractional quantum Hall ef- 
fects [10Hl5j . However, how to understand the real time 
evolution of an interacting system driven out of equilib- 
rium remains a great challenge in spite of intense efforts 
in recent years. The standard tool for the quantum field 
theory of nonequilibrium states is the Keldysh-Green's 
function technique. But the Keldysh techniques suffer 
from several shortcomings. First it is based on the Wick's 
theorem, and then the initial states must be limited to 
the equilibrium states of the quadratic Hamiltonian. Sec- 
ondly, the result from a perturbative expansion with re- 
spect to the interaction strength does not incorporate the 
physics of the problem when the Coulomb interaction is 
the largest energy scale. Finally, for some models, e.g., 
the Kondo model, the perturbation theory is plagued by 
the infrared divergence at low temperatures. 

For these reasons, new theoretical tools have been de- 
veloped to understand the interacting fermion systems in 
nonequilibrium. Numerical methods are developed like 
the real time Quantum Monte Carlo (16H19I. t he time- 
dependent numerical renormalization group [2CH22| . the 
scattering state numerical renormalization group [23j , the 
adaptive time-dependent density matrix renormalization 
group [2~4T - |26| and the non-equilibrium dynamical mean 



field theory [27j ■ The numerical method has the advan- 
tage that the result can be obtained in very high precision 
and the parameters of the model can be chosen arbitrar- 
ily. However, they also suffer the inconvenience that the 
large computer resources are required and the result does 
not incorporate the physics of the problem directly. 

The analytical methods are developed, which can sup- 
plement the shortcomings of the numerical methods, such 
as the scattering Bethe-ansatz (28[, the method based 
on integrability [29l |30| . the real-time renormalization 
group [3ll - [34| . the nonequilibrium flow equation (35r - 
|39(, and various approximation schemes building on the 
Green's function techniques (40ll4l|. However, no unique 
method is available which can cover all regimes of inter- 
est in different models. It is still necessary to study new 
methods. 

The Hamiltonian of a typical interacting fermion sys- 
tem can be expressed as 



k^ ,&2 



(1) 



where cj~ (ct) is the fermionic annihilation (creation) op- 
erator. And the normal ordering is with respect to the 
non-interacting ground state of the Fermi sea. In the 
traditional diagrammatic expansion the time evolution 
of a physical observable, denoted as (0(i)) where 0{t) 
satisfies the Heisenberg equation of motion 



i-0(t) = [0(t),H], 



(2) 
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is calculated by expanding the S-matrix and summing 
up a series of diagrams in a self-consistent way. In this 
paper we propose an alternate way to solve the Heisen- 
berg equation of motion by decomposing the observable 
operator O into the linear combination of the excitation 
operators which are defined as the operator A satisfying 
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the eigen equation 

[H,A] = XA. (3) 

This method circumvents the expansion of the S-matrix 
and then the contraction of field operators, and is appli- 
cable even the initial state is nontrivial. Furthermore, as 
we will show, this method is beyond the traditional per- 
turbation theory because the excitation energy A which 
appears in the exponent in the time function is modified 
by the interaction strength. 

In Sec. II we give the definition of the excitation oper- 
ators. In Sec. Ill we demonstrate how to construct the 
excitation operator and use it to calculate the time evo- 
lution of the creation operator in a toy model inspired by 
the Hubbard model. And the comparison between this 
model and the noncquilibrium flow equation approach 
which also contributes to solve the Heisenberg equation 
of motion is given. In Sec. IV we perform our method in 
single impurity Anderson model and calculate the evolu- 
tion of the current operator. Sec. V is a short summary. 

II. EXCITATION OPERATORS 

The excitation operator [42| of a Hamiltonian H is the 
operator A which satisfies 

[H,A]=xA, (4) 

where A is a real number denoting the excitation energy. 
We suppose that and \ip y ) are two arbitrary eigen- 
states of H with the eigen-energies E x and E y respec- 
tively. Then the excitation operator can be expressed 
as 

A x ,y = \lp x ){lpy\. (5) 

The corresponding excitation energy is X x>y — E x — E y . 
An observable operator O can be decomposed into the 
linear combination of the excitation operators, i.e., 

O = ^0 X ,yA x ,y, (6) 

where O x , y — {ip x \0\ip y ). According to Eq. 01 we have 

< '"' ■ \., ..,< '"' = e tX *^A x , y . (7) 

The solution of the Heisenberg equation of motion can 
be expressed as 

()ih X°'"'" v ( 8 ) 

x,y 

Once we know how to decompose an observable operator, 
we obtain the time evolution of it. 

In many-body problems it is generally difficult to find 
the eigen-states of the Hamiltonian and then calculate 



the excitation operator according to Eq. [5l An alternate 
way is to expand A into a series of products of field op- 
erators, substitute the expansion into Eq. 2] and decide 
the coefficients. For the Hamiltonian of Eq. [TJ we could 
express A as 

a = Mk£ l + Y Mk ' k : d l £k : 

+ m kk*i : 4i4s&i : + ■ • • > ( 9 ) 

where Mk, M^k and Myj^ki are t ne undetermined co- 
efficients. The commutator of the Hamiltonian with the 
products of TV field operators contains the products of 
(N + 2) operators, which indicate that A should be an 
infinite series. We need to truncate the series to get an 
approximated expression of A. 

Since we employ the perturbative expansion in the ex- 
pression of A, our method is also a kind of perturba- 
tive approach. While the direct perturbative approach 
for solving the Heisenberg equations of motion will gen- 
erate secular terms, which are uncontrolled at large 
timescales [3(1 [37| ■ By first decomposing the observable 
into the linear combination of the excitation operators, 
we avoid the secular divergences and make the results 
about the long-time behavior more explicit. 

III. A TOY MODEL HAMILTONIAN WITH 
DEGENERATE INTERACTIONS 

A. Excitation operators 

In this section, we use a toy model Hamiltonian 

H = Yl ek '- £ L & ka ■ Y '■ ^k'AaCkaCk'a ■ 

k,a k^k' ,<r 

+ 7T Y '■ B l'* 6 k* d >"9Cko ■ (10) 

to demonstrate how to construct excitation operators and 
use them to calculate the evolution of an observable op- 
erator. In Eq.[TU]the first term denotes the single particle 
energy and the second and third terms the interactions. 
Here a =f, i denotes the particle spin, and a the opposite 
to a. The normal ordering is with respect to the Fermi 
sea with chemical potential /i. The Fermi distribution 
function is denoted as — 6{n — efc). It is not difficult to 
find that this model Hamiltonian is inspired by the Hub- 
bard model. In fact it is the simplified Hubbard model in 
which most of the interaction terms are thrown off. Only 
the degenerate interaction term is left which is defined 
as the interaction keeping the sum of the single parti- 
cle energies of the two annihilation operators exactly as 
same as that of the two creation operators. For example, 
the degenerate interaction in the Hamiltonian of Eq. [1] 
should be the interaction terms Uk' 1 k' 2 k 1 k 2 '■ cjj' Ck 2 £k 1 ■ 
in which the relation €k x + £fc 2 = + £k' is satisfied. 
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We are interested in calculating e tHt cte lHt . We no- 
tice that cl increases the particle number of a state by 
one. Then it will be decomposed into the excitation oper- 
ators in which arbitrary term contains one more creation 
operator than the annihilation operator. We use the sym- 
bol A ka to represent the excitation operators. Again we 
could suppose that the first term of Aka is c fco . in the sim- 
plest case. The commutator of H with ct contains the 
terms : c\ a c\, s c v s : and : c\, a c\ & Ck> a ■■ The commuta- 



tor of H with : cl„ci,^c k ', 



and : c\, a c ks Ck' a '■ contains 



"ka^k 

one more term : c|»4' CT Cfc< a '■■ So we could suppose the 
expression of the excitation operator as 

Ak* = 4«r + 2 ( M fcV : AaA-s^'s ■■ 

k'^k 

+•••. (ii) 

In Eq. [TTJ we neglect the high order terms consisting of 
more than three field operators. 
The eigen-equation is written as 



[H, Aka] = ^kAka, 



(12) 



where Afe is the excitation energy. Substituting Eq. [TO] 
and [TT] in and comparing the coefficients before c kr7 be- 
tween the left and right sides of the equation, we get 

A fe = e k + U ]T (ML, + Ml k ,){-nl +n k ,). (13) 

k'^k 

But —n k , + rik' = whether rik> = or nk> = 1. Then 
the excitation energy becomes 

A fe = e k . (14) 

Comparing the coefficients before : c} ka $ k , s Ck*a ■ 

^k'tr^ka^k's '■ and : c ka c k/(7 c k >a : between the two sides 
of the eigen-equation for k ^ k' ', we get 

(A fe - e k )Ml k , =U (1 + Ml k ,(l - 2n k ,) 
+ M lk'( 1 ~n k - n k >)) , 

(Afe - £fe)M fe 2 fe, =U (1 + A4V(1 - n k - nfcO 
+Ml k ,(l-n k -n k >) 
+Ml k ,{n k -n k ,)) , 

(Afe - € k )Ml k , =UMl k ,(n k - n k ,). 

These equations are solved and the coefficients M are 
decided. When n k = ny we find that 



(15) 



1 



2n fc , 



(16) 



and Ml k , is an arbitrary number. What we want is a 
set of linearly-independent excitation operators. So we 
set M k k , — 0, since : c\ a c\, a c k i a : is itself an excitation 



operator as will be shown next. Eq. 1161 has two linearly- 
independent solutions: M k k , — 0, M kk , — 2n ; and 

^1 k' — 2n k ,-i > M kk , = 0. One could choose arbitrary 

one for Aka- We use an indicator function Xk,w = 0, 1 to 
denote our choice, and express the corresponding term in 
Aka as 



1 



2n k , - 1 
+ (1 - Xk,h 



Xk,k> 



st it 



C k'a C k'a '■ 



"ka^k'a 



^k'a^ks^k'a '■ 



(17) 



For n k 7^ n k j , the solution of Eq. [TS] is 

M x - 1 



Mi k > = o, 



Ml k , 



1 



(18) 



n k < - n k 

In summary the excitation operator can be expressed as 



Aka —C ka 



E 

n k ,=£n k 
1 



1 



ny - n k 



2n k , - 1 



^ka^k'a^k'a 



^ka^k'a^k'a ■ 



E 



i 



k'^k,n k , =n k 



2n k , - 1 



Xk,k' ■ cljl^Ck'a ■ 



+ (1-Xfe,fe') : A'aAa C k'S : ) 



(19) 



Different indicator functions will give different Aka- But 
we will show next that the choice of Xk,k' does not affect 
the last result of e lkt c{ a e~ liit . 

To decompose c ka we need to cancel the terms : 

Aa C l'a C k'a ■ c\, ■ J ks C k <a ■ and ! C^4'<A'<7 : in Aka. 

So we construct the excitation operators with them as 
the leading terms. To avoid ambiguity we use the sym- 
bol B to denote these operators. For arbitrary k ^ k' we 
suppose that 



B 



kk' 



Nl 



k.k 1 



st A t 



Cfc,=Cfe/ 



-N, 



k,k' 



^ka^k'a 



+ N lk> ■ c{, a ci-Ck>a ■ 



ka C k'a C k'a '■, 



(20) 



where we again neglect the high order terms consisting of 
more than three field operators. This supposition is rea- 
sonable because the commutator [H, B kk 'a\ does not con- 
tain the terms like single creation operator. The eigen- 
equation is now written as 



[H,B k k'a] - Af fe ,J5fcfc/ CT . 



(21) 



By comparing the coefficients between the left and right 
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TABLE I: The six types of B kk / a . operators 



Name 
B 1 
B 2 
B :i 
B 4 
B 5 
B 6 



Expression 

^ka^l'a^k'a '■ + '■ C fc , CT C fcff Cfc/ g 
C kcr C k' d C k'a '■ ~ ■ C k' o C ka C k' a 

■ C ka C k' a Ck ' " ' 

• C kcr C k'a C k'& ■ 

«t :t . , . at at 



k> cr C ka C k'a ■ + ■ C ka C k' a C k' a 



C k' o C ka C k' V 



k' <7 



Energy 
e k + 2U(l-2n k ,) 

e k + U(l - 2n k .) 
e k + U (n k - n k , ) 
e k -U(n k - ny ) 



sides of the equation we have 



i^k,k' 



(A 



B 

k,k' 



ek)Nl k , =U(Nl k ,(l-2n k ,) 

+Nl k ,(l-n k -n k ,)) , 

tk)Nl k , =U ({Nl k , + Nl k ,)(l - n k - nv) 
+Nh, k >(n k -n k >)) , 

ek)Nl k , =UNl k ,{n k -n k r). 



(22) 



In the matrix form of this equation we find that 
(N ktk , , N k k , , N k>k , ) is in fact the eigenvector of a (3 x 3) 
matrix. When n k = n k i Eq. [2U can be expressed in the 
matrix form as 



e k + U(l - 2n k > 
U{l-2n k .) 




Nlk> 
Nik' 



U(l-2n k ,) 
e k + U{l-2n k ,) 


Nl k , 



= A 



k,k' 



Nlk- 
Nik- 



(23) 



This matrix has three linearly- independent eigenvec- 
tors. Correspondingly we get three linearly-independent 
Bkk'a- When n k ^ n k , Eg. |2"2"1 becomes 



e k + U{l-2n k ,) 











Nl,k' 
Nik' 
Nik' 



A k,k' 



£fc U(n k -n k <) 
U(n k -n k >) e k 

Nik' 



Nl,k' 
Nik' 



(24) 



find that B\ k , a 



_ (4o-4'ff -~AsA'S c k's and the cor- 
responding excitation energy is e k + 2U. This operator 
annihilates an electron and simultaneously creats a spin 
singlet. This procedure costs an extra energy of 2U. 



B. The evolution of single creation operator 

The operator ct can be decomposed into the linear 
combination of the operators A klJ and B kk > a . The second 
order terms in A ka with n k — n k i can be canceled by 
Hk'cri B lk'a and pkk'a* and those with n k ^ n k> by 

B tk><r' B lk'a and ^fcfcV lt is eaS y t0 find 



C L -A ka - 



E 

k' ' ^k,n k i — rifc 



l 



2n k , - 1 V 2 



1 B 1 

~°kk'- 



+(Xk,k' - -j) B kk'a 



E 



2n k 



-B 



k k ' ■ 



2(n k , - rik) 



B 
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2(n k > - n k ) 



B 



kk' a 



(25) 



The expressions of e lHt A k<y e~ lHt and e lHt B kk /, 
got from Eq.[7|and the excitation energies listed in table I. 
Substituting back in the expression of A ka and B kk i a we 
get 



=e"**4«r + 



e*«n _ e 2»*t7(l-2n»/)) 



, E 2(2//,,. - I 

k'^k,n k f — rife 

x ( : 6 L4'&Cfe's = + : cl'As^k's ■) 
e 



E 



kt (l - e itU(l-2n k ,)j 

2nki _ I : c ka c k , 9 c k , s 



■ C k'o C kS Ck '° 



ie lttk sm[tU(n k , - n k )] _ ^ . t 
n k i - n k 

e it£k (l - cos[tt/(n fc - nfcQ]) . „ t at 
rife' - n fc 



c ko C k'o Ck 'o '■ 



(26) 



where we have used the identity n k +n k i = 1. This matrix 
contributes another three eigenvectors. Totally we obtain 
six types of B kk > a , which are all listed in table I with the 
corresponding excitation energies. 

In summary we get the excitation operators A ka and 
Bkk'a- We could also get the excitation operators of 
higher orders by similar analysis. Even we only keep 
the lowest order terms in the perturbative expansion 
of B kk i a , we get a nontrivial result. The interaction 
strength U enters the expression of the excitation energy 
of B kk i a . We conclude that B kk i a reflects the charac- 
teristics of the collective excitations. When n k i = we 



ifftat „-iHt 

L fc<T e 



is unique, inde- 



Note that the expression of e l 
pendent to the indicator function. 

In Eq. [26] the function of time is ?7-dependent, which 
indicates that our method involves the renormalization of 
the excitation energies and is beyond the traditional per- 
turbation theory. This equation is an operator identity. 
As the operator c ka (t) acting on e lHt \ip) where is an 

eigenstate, we get e iHt c ka \ip). So the physical meaning 
of Eq. [26] is that it shows how the system responses to 
a single particle excitation by producing the collective 
excitations. 
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It is worthwhile to mention that the nonequilibrium 
flow equation (3(| H3] is also a method designed for solv- 
ing the Heisenberg equation of motion by expanding the 
operator into a power series. There are three main dif- 
ference between our method and the flow equation ap- 
proach. First in the flow equation approach a unitary 
transformation is performed to the Hamiltonian, while 
in our method the Hamiltonian keeps invariant. Sec- 
ondly, we turns the problem into the diagonalization of 
the coefficient matrix instead of solving the differential 
equation in the flow equation approach. Finally we deal 
with the degenerate interaction in the same way as the 
non-degenerate ones, which suggests that this method 
is applicable in the model where the degenerate inter- 
action is overwhelming, e.g. in the Tomonaga-Luttinger 
model • While in the flow equation approach the de- 
generate interaction has a different scaling behavior in 
the flow. 



IV. SINGLE IMPURITY ANDERSON MODEL 

Next we consider the single impurity Anderson model, 
because it is important both in theory and experiment. 
Furthermore, many papers have contributed to study this 
model, so that we could compare our method with the 
others. We will use the excitation operators to calculate 
the time evolution of the current through the impurity. 

The Hamiltonian of the Anderson impurity model is 
expressed as 

kaa kaa 

~ dld a + Ud\d t d\d h (27) 

where a =\, \. denotes the electron spin and a = L,R 
the left and right leads. Here c\. aa and d\ are the single- 
electron creation operators in the leads and at the im- 
purity site respectively. At initial time the leads are in 
thermal equilibrium with chemical potentials re- 
spectively. The coupling between leads and the impurity 
is switched on at time t = 0. 

In the pre-diagonalization basis [43|, HiJ , the Hamilto- 
nian can be expressed as 

ka sij 

+ UB S[S , 2S1S2 : cl, it cl, 2i c S2l c Slt :, (28) 

s' 1 s' 2 s 1 s 2 

where B s ' iS ' slS2 is the abbreviation of B s ' i B Sl B s > 2 B S2 
and the normal ordering is with respect to the ini- 
tial state. The (anti) symmetric operator is defined as 
Ck± a — -^(ckLa ± ctRa) and the hybridization operator 

as c sa = ^2 k j-^-B s 6k+a + B s d a . Here the coefficient 
B s = . v and the linewidth T = pnV 2 , where p is 



the density of states of the lead. In the Anderson impu- 
rity model we will neglect the degenerate interaction by 
ordering B s ^ SlS2 = as e s ' t + e s > 2 = e Sl + e S2 . This is 
reasonable in the Anderson impurity model, since the co- 
efficient B s — > in the thermodynamic limit. The Fermi 
function in hybridization basis has sharp edges in the 
thermodynamic limit, i.e., 

(c; v c SCT )o = S s , s 'n s , (29) 

where n s = ±(0(^ - e s ) + 9(~^f - e,)). 

We are interested in the evolution of the current oper- 
ator 

= 1 dN Lt dN m 
t 2 [ dt dt ' 

^Jfl^t + k., (30) 

where the antisymmetric field operator £fc_f is just the 
excitation operator of the Hamiltonian and its evolution 
can be written as e lHt Ck-fe~ %m = e~' l€kt £k-^- So the 
point is how to calculate e lHt c f s ^e~ lHt , which needs us to 
decompose into the excitation operators. 

Because in the Hamiltonian (|2"T|) the total electron 
number and spin are conserved, we suppose the excita- 
tion operator as 

a = J2 M ^U+ E K4°* :d lA d ^ : 

+ E M s[s' 2 s 2 ■■ ^ it ^At : • (31) 

Here we truncate the series of A and neglect the higher 
order terms consisting of more than three field operators. 
In this approximation we can obtain the current to order 
U 2 . 

Substituting Eq. ED into Eq. @] and comparing the co- 
efficients between the terms at the left and right sides of 
the equation, we get 

e s Ml + U ]T B ss , iS , 2S2 Q s{s , 2S2 M\ s , 2S2 =\Ml 
(e s; + e s , 2 - e S2 )M 2 , s , S2 + r ^ /i... I/; 

S 

+ U ^2 B s' 1 t lS2 t 2 {n tl - n t2 )M 2 iS , 2t2 
t±t 2 

+Uj2 B tit 2S ' lS ' 2 {l - n tl - n t2 )Ml t2S2 
t±t 2 

+U B tlt2S2S2 (n tl - n t2 )Ml s , it2 

tlt 2 

+UY / B tl t 2 s 2 s 2 (n t2 - n tl )M\ tit2 = AM S 2 , S , S2 , 
tits 

+£s> 2 ~ t S2 )M 3 s , iS , 2S2 + U^B s , 2S2tlt2 {n tl - n ta ) 

tits 

xM 2 , t2ti =\M3, s , S2 , (32) 
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where Q s > lS < 2 s 2 = n s tn s ' 2 - n s <n S2 - n s > 2 n S2 + n S2 . We 
express above equations in a matrix form 

HM = \M, (33) 

where M = (M},M 2 , , , M 3 , , ) is a vector. The di- 

mension of the vector M is (A + 27V 3 ), and it is the 
eigenvector of a (N + 2N 3 ) x (TV + 27V 3 ) matrix H, where 
A is the number of possible values of s, i.e., the number 
of levels in each lead. And A is the corresponding eigen- 
value. To find {Ml , M 2 s , S2 , M 3 s , S2 ) satisfying Eq. [32] is 
equivalent to diagonalize the matrix %. Now the prob- 
lem changes into the problem of diagonalizing a matrix 
whose dimension is much less than the dimension of the 
original Hamiltonian. The diagonal elements of T-L are 
{/-independent, while the non-diagonal elements are in 
order U . When U is small, we could employ perturbation 
theory to calculate the eigenvalues and the eigenvectors 
of H. By setting U = 0, we get the eigenvectors and 
the corresponding eigenvalues in zeroth order. There are 
three types of eigenvectors. They are 

Ml ~ {Ml = 1, M\ s , S2 = 0, M\ a , 2S2 = 0) (34) 
with the eigenvalue e s and 



Mi s , S2 ~ {Ml=0,M t 



1,M5 



1 "'S S'So 



0), 



M\ S , 2S2 ~ {Ml = 0,M 2 , S , S2 = 0,M 3 , 4S2 = 1) (35) 

with the eigenvalue + e s ^ — e S2 . For simplicity, 
we express the components of the eigenvector Ail as 
= <y> <%, S2 = M\%, S2 = 0, where the 
sup(sub)script before comma is the index of eigenvec- 
tors, and that after comma the index of components in 
each eigenvector. The same scheme is applied to the vec- 
tors M 1, „/ „ and M 3 , „, „ . For all of the three types of 

s l s 2 2 S 1 S 2 S2 A 

M the first and second order corrections to the eigenval- 
ues will go to zero in the thermodynamic limit when the 
level spacing and the coefficient B s go to zero. And the 
corrected eigenvectors are calculated to be 



U 2 B S B S ,T{D) 
e s - e s >){D - e a 



0(t/ 3 ), 



M 1 ' 2 , . 

S,S',S' S-2 



0{U% 



+ 0{U 2 ), (36) 



where D = ey + e s > — e s , is the eigenvalue of M 2 , , 

12 ^1^2* 



and Mi .> and 



T{D) = Y,Q^(^ 2 -o)BlB 2 s 



Corresponding to the three types of M , the three types 
of excitation operators can be expressed as 



u = 4 t + E <14> 



5^S 2 S2 



A 2 

■ rL s' 1 S 2 S2 



l s' 1 S 2 S 2 ,t' 1 t 2 t 2 ' tit * 2 4- * 2j - 

V A/f 2 ' 3 • r f r f 7v * • 

l s' lS2 s 2 ,tit 2 i 2 ■ L tit t 2 t * 2 T ■' 

i 1 ^ 2 ^ 2 



and 



A 3 , , 

s[s' 2 s 2 



= V M 3 ' 1 c f 

t 1 t 2 i 2 

+ 1 6 sit^S 2 t^2t 1 



t 2 t 2 ' ^t'jt^jA 2 ^ 



V A4 3 ; 3 , , , 

/ 4 s 1 s 2 S2,i 1 i 2 
i / 1 t 2 i2^s^s 2 S2 



(38) 



In the matrix form these equations can be expressed as 



A l lS ' 2 s 2 I = M 



: c *it c tu c ' 2j - 
V : S tit S * 2 T^ t2t 



(39) 



where M is a (A + 2A 3 ) x (A + 2A 3 ) matrix with the 
elements M',, M?) 3 , , etc.. We solve this system 

s.s 1 s 1 s 2 s 2 .t i t 2 t 2 1 J 

of linear equations and get the expression of cj^, i.e., 

+ E (^ _1 )iS^^*- ( 4 °) 

s^s 2 s 2 

We notice that the elements of the matrix {M — 1) are 
all in order U , where 1 is the identity matrix. Then for 
small U the inverse of M can be expressed by a Neumann 
series 



xB 2 {e s[ +e s , 2 -D). 



(37) 



Here we only list the components of M which will be 
used in the calculation of the current to second order of 
U. 



M' 1 = 1 - {M - 1) + {M - l) 2 - 



(41) 



The excitation energy of Al is e s , and that of A 2 ,^ s ,^ S2 
and A 3 , , . is D = e a < + e s < — e S9 . Then the evolution of 

S,S 2 S 2 S l 6 2 * 2 
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the hybridization operator can be expressed as 

s" s' 

+ Y CM- 1 ) 1 '' > M 2 ; 1 , „e lDt 

+ V (.M" 1 ) 1 ' 3 , , M 3 ; 1 , „e lDt )c\,, + --- . 

/ j \ < / S,S^S 2 S2 S^S 2 S2,S" / S"j 

S f 2 S2 

(42) 

In this expansion we only keep the lowest order term 
with single creation operator and neglect the higher or- 
der terms, because the higher order terms are normal 
ordered and will not contribute to the expectation value 
of the current operator. Substituting Eq. I4T1 and l36l into 
Eq.[42j we will obtain the coefficient before $ a „^ to second 

order of U. After getting the expression of e lHt cl^e~ lHt , 
we could directly express the evolution of the current op- 
erator according to Eq. [3D] as 

e im^ e -iHt = ^ 7s (t) e -^*4 t c fe _ t + h.c. + ■■■ , (43) 

sk 

where 

7 , (t) . ^, + 'ZMl EnD)B i 

si,D 

X _(e s - D)(e Sl - D) + (e s - e Sl )(e Sl - D) ' 

(44) 



Eq. [44] reproduces the result in Ref. [3 EH , which in- 
dicates that our method is a well-defined perturbative 
method in this context. However, we would like to men- 
tion that by employing the numerical routines as diago- 
nalizing the matrix % our method has the possibility to 
go beyond the perturbation theory. 



V. CONCLUSIONS 



We develop a method for calculating the real time evo- 
lution of an observable operator satisfying the Heisen- 
berg equations of motion. The point of this method is 
to construct the excitation operators of the Hamiltonian 
and then decompose the observable operator into the lin- 
ear combination of the excitation operators. We use this 
method to calculate the evolution of single creation oper- 
ator in a toy model Hamiltonian inspired by the Hubbard 
model and the evolution of the current operator in sin- 
gle impurity Anderson model. We expect to perform our 
method in a lattice system, e.g. the Hubbard model, in 
future studies. 
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